Import data

languages <- read_csv("voicing-effect/stimuli/languages.csv")
## Parsed with column specification:
## cols(
##   speaker = col_character(),
##   language = col_character()
## )
words <- read_csv("voicing-effect/stimuli/nonce.csv")
## Parsed with column specification:
## cols(
##   item = col_integer(),
##   word = col_character(),
##   ipa = col_character(),
##   c1 = col_character(),
##   c1phonation = col_character(),
##   vowel = col_character(),
##   anteropost = col_character(),
##   height = col_character(),
##   c2 = col_character(),
##   c2phonation = col_character(),
##   c2place = col_character(),
##   language = col_character()
## )
columns <- c(
    "speaker",
    "seconds",
    "rec.date",
    "prompt",
    "label",
    "TT.displacement.sm",
    "TT.velocity",
    "TT.velocity.abs",
    "TD.displacement.sm",
    "TD.velocity",
    "TD.velocity.abs"
)

aaa.files <- list.files(
    path = "./voicing-effect/results/ultrasound",
    pattern = "*-tongue-cart.tsv",
    full.names = TRUE
)

tongues <- read_aaa(
    aaa.files,
    columns, 
    na.rm = TRUE
) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor) %>%
    group_by(speaker) %>%
    mutate(
        X.re = rescale(X),
        Y.re = rescale(Y)
    ) %>%
    ungroup() %>%
    mutate(
        vowel.ord = ordered(vowel, levels = c("a", "o")),
        c2place.ord = ordered(c2place, levels = c("coronal", "velar")),
        c2phonation.ord = ordered(c2phonation, levels = c("voiceless", "voiced"))
    )
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_1 = col_character(),
##   Y_1 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Joining, by = "speaker"
## Joining, by = c("word", "language")
max <- tongues %>%
    filter(label %in% c("max_TT", "max_TD"), vowel != "u") %>%
    arrange(rec.date, fan.line) %>%
    create_event_start("rec.date")

max_it_12 <- max %>%
    filter(speaker %in% c("it01", "it02"))

max_pl_23 <- max %>%
    filter(speaker %in% c("pl02", "pl03"))

max_pl_34 <- max %>%
    filter(speaker %in% c("pl03", "pl04"))

Exploration

max %>%
    filter(c2place == "coronal", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

max %>%
    filter(c2place == "velar", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

max %>%
    filter(c2place == "coronal", language == "polish") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

max %>%
    filter(c2place == "coronal", language == "polish") %>%
    ggplot(aes(X.re, Y.re, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

max %>%
    filter(c2place == "velar", language == "polish") %>%
    ggplot(aes(X.re, Y.re, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

Testing

Comparing fixed effects models works best with ML, otherwise you can use fREML with AIC, but if there is an AR1 model cannot use AIC.

Polish

pl.m1 <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl.m1)

pl.m1.ar <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML",
    rho = rho,
    AR.start = max_pl_23$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X.re, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.2474     0.1778   1.392   0.1641  
## X.re          0.6437     0.3900   1.651   0.0989 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F  p-value    
## s(X.re)                         6.804   7.090 32.225  < 2e-16 ***
## s(X.re):c2phonation.ordvoiced   7.958   8.606  7.478 1.56e-08 ***
## s(X.re):c2place.ordvelar        8.825   8.965 85.509  < 2e-16 ***
## s(X.re):vowel.ordo              8.717   8.935 19.070  < 2e-16 ***
## s(X.re,rec.date)              481.302 540.000 19.552  < 2e-16 ***
## s(X.re,speaker)                 5.570   8.000 32.388  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 592/593
## R-sq.(adj) =  0.958   Deviance explained = 96.3%
## fREML = -6475.4  Scale est. = 0.0015912  n = 4212
pl.m1.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl.m1, pl.m1.null)
## pl.m1: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X.re, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## pl.m1.null: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2place.ord, 
##     bs = "cr") + s(X.re, by = vowel.ord, bs = "cr") + s(X.re, 
##     rec.date, bs = "fs", xt = "cr", m = 1, k = 5) + s(X.re, speaker, 
##     bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Chi-square test of fREML scores
## -----
##        Model     Score Edf  Chisq    Df   p.value Sig.
## 1 pl.m1.null -6456.357  14                            
## 2      pl.m1 -6475.358  16 19.001 2.000 5.596e-09  ***
## 
## AIC difference: -20.21, model pl.m1 has lower AIC.
plot_smooth(
    pl.m1,
    view = "X.re",
    plot_all= "c2phonation.ord",
    cond = list("c2place.ord" = "coronal"),
    rug = FALSE
)
## Summary:
##  * X.re : numeric predictor; with 30 values ranging from 0.001193 to 1.000000. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14. 
##  * speaker : factor; set to the value(s): pl02.

plot_diff(
    pl.m1,
    view = "X.re",
    comp = list(c2phonation.ord = c("voiceless", "voiced")),
    cond = list("c2place.ord" = "coronal")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.001193 to 1.000000. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14. 
##  * speaker : factor; set to the value(s): pl02.

## 
## X.re window(s) of significant difference(s):
##  0.001193 - 0.031460
##  0.243328 - 0.283684
##  0.868844 - 1.000000
plot_gamsd(
    pl.m1,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.001193 to 1.000000. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14. 
##  * speaker : factor; set to the value(s): pl02.

plot_gamsd(
    pl.m1.ar,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.001193 to 1.000000. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14. 
##  * speaker : factor; set to the value(s): pl02.

PL02

pl02_max <- filter(max, speaker == "pl02")

pl02.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl02.m1)

pl02.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "fREML",
    rho = rho,
    AR.start = pl02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl02.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.59019    0.60248  -4.299 1.79e-05 ***
## X            0.29360    0.03819   7.688 2.23e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df       F  p-value    
## s(X)                         7.555   7.894  69.913  < 2e-16 ***
## s(X):c2phonation.ordvoiced   8.011   8.702   7.806 1.95e-09 ***
## s(X):c2place.ordvelar        8.883   8.986 316.209  < 2e-16 ***
## s(X):vowel.ordo              8.201   8.788  20.088  < 2e-16 ***
## s(X,rec.date)              189.728 300.000   9.633  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 342/343
## R-sq.(adj) =   0.95   Deviance explained = 95.4%
## fREML =   5616  Scale est. = 4.5043    n = 2419
pl02.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl02.m1, pl02.m1.null)
## pl02.m1: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl02.m1.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of fREML scores
## -----
##          Model    Score Edf  Chisq    Df   p.value Sig.
## 1 pl02.m1.null 5638.118  11                            
## 2      pl02.m1 5615.978  13 22.140 2.000 2.425e-10  ***
## 
## AIC difference: -37.43, model pl02.m1 has lower AIC.
plot_gamsd(
    pl02.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -26.245000 to 59.219800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14.

PL03

pl03_max <- filter(max, speaker == "pl03")

pl03.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl03.m1)

pl03.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML",
    rho = rho,
    AR.start = pl03_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl03.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML",
    rho = rho,
    AR.start = pl03_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl03.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.15251    0.84915  -6.068 1.63e-09 ***
## X            0.42071    0.05257   8.004 2.35e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.433   7.793 52.044  < 2e-16 ***
## s(X):c2phonation.ordvoiced   6.740   7.654  4.817 1.54e-05 ***
## s(X):c2place.ordvelar        8.755   8.938 33.136  < 2e-16 ***
## s(X):vowel.ordo              8.770   8.941 20.139  < 2e-16 ***
## s(X,rec.date)              209.072 235.000 26.599  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 277/278
## R-sq.(adj) =  0.969   Deviance explained = 97.3%
## fREML = 3995.5  Scale est. = 2.8743    n = 1793
pl03.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl03.m1.ar, pl03.m1.ar.null)
## pl03.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl03.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl03.m1.ar.null 3335.866  11                         
## 2      pl03.m1.ar 3330.848  13 5.018 2.000   0.007  **
## Warning in compareML(pl03.m1.ar, pl03.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.611727, rho2 = 0.611727).
plot_gamsd(
    pl03.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

PL04

filter(max, speaker == "pl04") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ c2place)

filter(max, speaker == "pl04", X > -30) %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ c2place)

pl04_max <- filter(max, speaker == "pl04", X > -30)

pl04.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04.m1)

pl04.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "ML",
    rho = rho,
    AR.start = pl04_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl04.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "ML",
    rho = rho,
    AR.start = pl04_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.28296    0.38386   8.553   <2e-16 ***
## X           -0.03567    0.03469  -1.028    0.304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.538   7.847 31.856  < 2e-16 ***
## s(X):c2phonation.ordvoiced   1.006   1.009  0.296    0.589    
## s(X):c2place.ordvelar        8.539   8.876 26.159  < 2e-16 ***
## s(X):vowel.ordo              7.975   8.620 10.694 7.37e-15 ***
## s(X,rec.date)              197.731 230.000 25.617  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 272/273
## R-sq.(adj) =  0.975   Deviance explained =   98%
## -ML = 1757.9  Scale est. = 0.72325   n = 1228
pl04.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
# compareML(pl04.m1, pl04.m1.null)
compareML(pl04.m1.ar, pl04.m1.ar.null)
## pl04.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl04.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04.m1.ar.null 1758.008  11                         
## 2      pl04.m1.ar 1757.861  13 0.148 2.000   0.863
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.518140, rho2 = 0.518140).
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): Only small difference in ML...
plot_gamsd(
    pl04.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -29.959700 to 32.556800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:58:11.

acf_plot(resid(pl04.m1.ar), split_by=list(pl04_max$rec.date))

acf_resid(pl04.m1.ar, split_pred = "AR.start")

pl04_a <- filter(max, speaker == "pl04", vowel == "a")

pl04.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04.m1)

pl04.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "ML",
    rho = rho,
    AR.start = pl04_a$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl04.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "ML",
    rho = rho,
    AR.start = pl04_a$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.83223    0.46029   1.808   0.0711 .
## X            0.02144    0.04680   0.458   0.6470  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df      F p-value    
## s(X)                        7.479   7.786 24.294  <2e-16 ***
## s(X):c2phonation.ordvoiced  3.384   4.190  0.553   0.575    
## s(X):c2place.ordvelar       8.019   8.623 12.720  <2e-16 ***
## s(X,rec.date)              95.926 116.000 29.624  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 148/149
## R-sq.(adj) =  0.987   Deviance explained = 98.9%
## -ML =   1101  Scale est. = 0.95795   n = 689
compareML(pl04.m1.ar, pl04.m1.ar.null)
## pl04.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## pl04.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04.m1.ar.null 1101.707   9                         
## 2      pl04.m1.ar 1101.011  11 0.696 2.000   0.499
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.455701, rho2 = 0.455701).
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): Only small difference in ML...

If only /a/ is used, no significance at all. With /o/ and X > -30, no significance, but with all X some significant difference at the very left end of the tongue.

plot_gamsd(
    pl04.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -41.892500 to 32.556800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:59:49.

PL03-PL04

pl34.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl34.m1)

pl34.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl34.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl34.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5) + s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -4.8776     3.6453  -1.338    0.181
## X             0.2787     0.2154   1.294    0.196
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         6.853   7.055 52.880  < 2e-16 ***
## s(X):c2phonation.ordvoiced   5.106   6.157  1.935   0.0706 .  
## s(X):c2place.ordvelar        8.572   8.892 38.777  < 2e-16 ***
## s(X):vowel.ordo              7.293   8.174  7.801 1.55e-10 ***
## s(X,rec.date)              354.481 469.000 10.193  < 2e-16 ***
## s(X,speaker)                 5.618   8.000 34.261  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 522/523
## R-sq.(adj) =  0.959   Deviance explained = 96.4%
## -ML = 5944.5  Scale est. = 2.5902    n = 3130
compareML(pl34.m1.ar, pl34.m1.ar.null)
## pl34.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5) + s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## pl34.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl34.m1.ar.null 5947.104  14                         
## 2      pl34.m1.ar 5944.513  16 2.591 2.000   0.075
## Warning in compareML(pl34.m1.ar, pl34.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.609938, rho2 = 0.609938).
## Warning in compareML(pl34.m1.ar, pl34.m1.ar.null): Only small difference in ML...
plot_gamsd(
    pl34.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar", speaker = "pl04")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -43.836000 to 70.078500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40. 
##  * speaker : factor; set to the value(s): pl04.

The models with X, Y are not reliable because the tongue size differs a lot between speakers.

pl34.m1.re <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl34.m1.re)

pl34.m1.ar.re <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl34.m1.ar.re.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl34.m1.ar.re)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.07931    0.03880   2.044    0.041 *  
## X.re         1.10388    0.08357  13.209   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F  p-value    
## s(X.re)                         7.433   7.784 38.617  < 2e-16 ***
## s(X.re):c2phonation.ordvoiced   5.673   6.770  2.625   0.0127 *  
## s(X.re):c2place.ordvelar        8.417   8.820 14.290  < 2e-16 ***
## s(X.re):vowel.ordo              8.418   8.822 11.068 6.22e-16 ***
## s(X.re,rec.date)              408.292 469.000 12.759  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 512/513
## R-sq.(adj) =  0.964   Deviance explained = 96.9%
## -ML = -6142.8  Scale est. = 0.0010642  n = 3130
compareML(pl34.m1.ar.re, pl34.m1.ar.re.null)
## pl34.m1.ar.re: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## pl34.m1.ar.re.null: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2place.ord, 
##     bs = "cr") + s(X.re, by = vowel.ord, bs = "cr") + s(X.re, 
##     rec.date, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                Model     Score Edf Chisq    Df p.value Sig.
## 1 pl34.m1.ar.re.null -6138.175  11                         
## 2      pl34.m1.ar.re -6142.819  13 4.643 2.000   0.010  **
## Warning in compareML(pl34.m1.ar.re, pl34.m1.ar.re.null): AIC is not
## reliable, because an AR1 model is included (rho1 = 0.590508, rho2 =
## 0.590508).
## Warning in compareML(pl34.m1.ar.re, pl34.m1.ar.re.null): Only small difference in ML...
plot_gamsd(
    pl34.m1.ar.re,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.015449 to 0.982036. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

plot_gamsd(
    pl34.m1.ar.re,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.015449 to 0.982036. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

Italian

it.m1 <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it.m1)

it.m1.ar <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML",
    rho = rho,
    AR.start = max_it_12$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X.re, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.25090    0.06259   4.008 6.26e-05 ***
## X.re         0.80000    0.21532   3.715 0.000207 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df       F  p-value    
## s(X.re)                         5.516   6.082   3.409  0.00394 ** 
## s(X.re):c2phonation.ordvoiced   5.106   6.209   5.491 1.33e-05 ***
## s(X.re):c2place.ordvelar        8.706   8.942 101.743  < 2e-16 ***
## s(X.re):vowel.ordo              5.286   6.405   8.554 1.45e-09 ***
## s(X.re,rec.date)              177.731 415.000   1.097  < 2e-16 ***
## s(X.re,speaker)                 6.607   8.000  74.894  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 467/468
## R-sq.(adj) =  0.941   Deviance explained = 94.5%
## -ML = -6713.9  Scale est. = 0.0017165  n = 3165
it.m1.ar.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
#        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML",
    rho = rho,
    AR.start = max_it_12$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
it.m1.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it.m1.ar, it.m1.ar.null)
## it.m1.ar: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X.re, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## it.m1.ar.null: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2place.ord, 
##     bs = "cr") + s(X.re, by = vowel.ord, bs = "cr") + s(X.re, 
##     rec.date, bs = "fs", xt = "cr", m = 1, k = 5) + s(X.re, speaker, 
##     bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##           Model     Score Edf  Chisq    Df   p.value Sig.
## 1 it.m1.ar.null -6702.103  14                            
## 2      it.m1.ar -6713.863  16 11.761 2.000 7.805e-06  ***
## Warning in compareML(it.m1.ar, it.m1.ar.null): AIC is not reliable, because
## an AR1 model is included (rho1 = 0.762798, rho2 = 0.762798).
plot_gamsd(
    it.m1.ar,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.000000 to 0.995850. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52. 
##  * speaker : factor; set to the value(s): it01.

IT01

it01_max <- filter(max, speaker == "it01")

it01.m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it01.m1)

it01.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "ML",
    rho = rho,
    AR.start = it01_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it01.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.73141    0.30682  -15.42   <2e-16 ***
## X            0.76532    0.02666   28.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df       F  p-value    
## s(X)                        7.764   7.944 181.554  < 2e-16 ***
## s(X):c2phonation.ordvoiced  5.037   6.136   9.480 2.06e-10 ***
## s(X):c2place.ordvelar       8.822   8.979 183.531  < 2e-16 ***
## s(X):vowel.ordo             6.862   7.808  11.566 1.19e-15 ***
## s(X,rec.date)              92.680 225.000   1.976  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 267/268
## R-sq.(adj) =  0.979   Deviance explained =   98%
## -ML = 3488.2  Scale est. = 4.0701    n = 1932
it01.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "ML",
    rho = rho,
    AR.start = it01_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it01.m1.ar, it01.m1.ar.null)
## it01.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it01.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it01.m1.ar.null 3509.547  11                            
## 2      it01.m1.ar 3488.231  13 21.316 2.000 5.527e-10  ***
## Warning in compareML(it01.m1.ar, it01.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.745137, rho2 = 0.745137).
plot_gamsd(
    it01.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -40.731200 to 54.425500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52.

IT02

it02_max <- filter(max, speaker == "it02")

it02.m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it02.m1)

it02.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "ML",
    rho = rho,
    AR.start = it02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it02.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.67360    0.52437  -3.192  0.00145 ** 
## X            0.65159    0.03157  20.640  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df      F  p-value    
## s(X)                        6.386   7.189 57.403  < 2e-16 ***
## s(X):c2phonation.ordvoiced  6.589   7.613  4.927 9.93e-06 ***
## s(X):c2place.ordvelar       8.630   8.910 52.382  < 2e-16 ***
## s(X):vowel.ordo             5.937   7.020  9.062 6.16e-11 ***
## s(X,rec.date)              72.997 185.000  1.401  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 227/228
## R-sq.(adj) =  0.916   Deviance explained = 92.2%
## -ML = 2803.6  Scale est. = 9.778     n = 1233
it02.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "ML",
    rho = rho,
    AR.start = it02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it02.m1.ar, it02.m1.ar.null)
## it02.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it02.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df   p.value Sig.
## 1 it02.m1.ar.null 2813.488  11                           
## 2      it02.m1.ar 2803.592  13 9.897 2.000 5.034e-05  ***
## Warning in compareML(it02.m1.ar, it02.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.735050, rho2 = 0.735050).
plot_gamsd(
    it02.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -44.949000 to 62.666500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:45:14.